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A KALMAN FILTER FOR A POISSON SERIES WITH COVARIATES 

AND LAPLACE-APPROXIMATION INTEGRATION 

D. P. Gaver 
P. A. Jacobs 



0. INTRODUCTION 

The Poisson model is an initial idealized, but plausible, off-the-shelf tool 
for representing point-process data of nearly any kind, cf. Feller (1966) and 
Cox and Lewis (1968). However, to be more descriptive, and even predictive, 
representation of non-homogeneity in space or time may be needed. For 
instance the occurrence of rare events in space, such as the occurrence of 
extreme heights, i.e. above a level of Arctic ice along a transsect or encounters 
with ice keels along a submarine track at constant (deep) depths, may well 
appear roughly Poisson. A better description of these events may require 
more detail than a simple mean or rate: some account of regional and 
seasonal variation could be needed for true accuracy; for instance the 
intervention of natural gaps along a path in Arctic ice will occur if the latter 
crosses leads (open water amidst an ice pack). For another example, demands 
for spare parts in a logistics system often are roughly Poisson, or compound 
Poisson, but with a mean or rate that changes erratically but slowly in time. 
Demands for communication and computer facilities exhibit a similar 
temporal pattern, and there are a great many other examples. To summarize, 
variation in a fundamental Poisson rate or mean is very often encountered in 
practice; it is possible that if this variation is slow-moving or persistent 
enough it can be exploited for short-term forecasting. 
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In this paper we study a model-based procedure for forecasting in the 
kind of environments described. The model introduced allows the mean or 
rate of the Poisson process to be itself a random process; the exponential of an 
AR/1 autoregressive process. In addition the rate is influenced by a 
covariate. We then recursively update the parameter estimates using an 
approximation based on the Laplace method ; cf. de Bruijn (1958). The approach 
resembles that of Delampady, Yee and Zidek (1991), but frankly heuristic 
methods are used to estimate certain of the underlying parameters. The 
methodology is checked against simulated data with encouraging results. 

1. FORMULATION 

Consider the following Poissonian model for count data: 

f 

P{Y, =y,|p,,i,,P,Vyci.-.yi- 1 } - exp{-Jr,e I,p * |1 < ) (1.1) 

iff 

where • • 

(1-2) 

with {co f } independent normal/Gaussian random variables with mean 0 and 
variance W t independent of {\i;; i < f-1), {Y,-; i < t-1}, {x t } r {h t } and p. Another 
hierarchical time series model for count data can be found in Harvey et al. 
(1989). 

The purpose of this paper is to suggest a Kalman filter-like procedure to 
produce successive estimates of P and |i t as new data becomes available. The 
procedure is based on a Laplace approximation to an integral. A Bayesian 
approach to a similar problem is being investigated by Delampaday et al. 
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(1991). A Bayesian approach to time series can be found in West et al. (1989) 
and for a more recent computational approach in Carlin et al. (1991). 



2. AN APPROXIMATE UPDATING PROCEDURE 

Assume that the posterior distribution of (JJ, ji (1 ) given {y,-, i< t- 1} is 
bivariate normal with mean ( b t _ v m M ), Var[|J] = x t _ v Var[p. (1 ] = C t _j and 
Corr[p,n M ] = p M . 

Since it is known that 

Rt =aji,.i +<o f . 



the prior distribution of (fJ,Ht) is bivariate normal with mean (£> f _ 1 , aw H ), 
Var[p] = t m , Var[n,] = R t = aC t _ x + W, and Corr[(l,n,] = r t = ap t _ x VC/_i/R/. 

The forecast/prediction distribution of in terms of data up to t-1 and 
the covariate value at t is 



P {Y ( = y f |Y, = y 2 -,z < t - l,x if i < f,fy , i < f} 

, '[ e' x ' b '\Y‘ 



exp 



iMf 1 



y t '- 2^T ( . 1 R f (l - r f 2 ) 

(b-bf) 2r t (b - bf)(z ■ mf) (z -mf) 






Rt 



■dzdb. (Zl) 



where bf = b t . i and mf = am t .j. 



We now approximate the integral by the Laplace method; cf. Easton 
(1991), Cox and Hinkley (1974), de Bruijn (1958). Let the exponent of the 
integrand be 
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( 2 . 2 ) 



g(b,z) = -e^^ht + y t [x t b + z] + y t \nh t - |(l- r f 2 ) 1 



Vi ,/WV * n " R, 



+ K 



where K is a constant. 



#-*(6,2) - y, ■ 



dz 



1 * ri 



( z : m t ) r . ( b ' b f ) 

Rt 1 >/ T f -l 






^^(^ 2 ) = [y t ■ e Xtb+z h t )x t - [l - r}\ 



-1 



( b - b t) r 4 2 • 

T f-i V T f-i V^7 



&■ 



-$(&,z) = -e( Xtb+z) h t x} -[l -r, 2 ] 



•1 1 



h-\ 



dbdz 



g(b,z)= -e^ Xtb+z) h t x t +(r f )(l - r f 2 ) V^R,) 
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Use a Newton procedure to solve the system of equations 



o-A, 

db 6 



for mt and bt . Solve for x lr C t and p t using the relation 



(2.3) 

(2.4) 

(2.5) 

( 2 . 6 ) 

(2.7) 

(2.8) 
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Pt-JC t T t 




r e 2 

— “T & 




a 2 

_ g 










a b 2 


mt,b t 


dbdz 


m t ,b t 


1 


C t 


| 


a 2 

m dbdz ^ 


m tM 


a 2 

2 % 
dz L 


mtM „ 



The posterior distribution for (P41O given yt, D t ~\ is approximated by a 
bivariate normal distribution with mean ( bt , mt) and variance-covariance 
matrix 




( 2 . 10 ) 



2.1 Summary of the Newton Procedure 

0 . Start with estimates of the parameters of the prior bivariate distribution 

of that is, estimates of the mean (b t -\, am t - 1), Corr(P,p. f ) = r t , Var(P) = r t -\ 

Var[p. f ] = R f . Set b° t = b^_ v m° t = am t _ v 

1 . Solve the system of linear equations (in b and m ) 




for b, m . 

2. Ifmax^|b ( ° -b |/b),(|m® -m|/m)J< 0.001, set m t = m and b t = b and go to 

3 . Otherwise set b° t = b and m ° t = m and return to Step 1 unless Step 1 has been 
returned to 49 times in which case set m t = m° t and b t = bt ; go to 3 . 

3 . Return ( bt,ntt ) as the estimate of the posterior mean of the bivariate 
normal distribution of (P,p f )- 
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4. To find the variance-covariance matrix of the posterior distribution of 
(p,|X f ) evaluate 



db‘ 



dzdb 



t ) 






dzdb 
=>2 



8i b t>™t) 



“2 Stow) 

dz 



set it equal to 

PtJhJCt 

S>t4hjc"t Q 

and solve for x t , p t and Q. 

2.2 Summary of Kalman Procedure 

In summary the approximate Kalman procedure is as follows 

0. Start with the parameters of the (approximate) posterior bivariate 
normal distribution of having mean (bt- Corr(p,p.f_i) = p t -\, 



Var[p] = tm Var[p,_i] - C M - 



1. Update Corr(p,p.<_i) and Var[^i(.i] using (1.2) to obtain the prior 
bivariate normal distribution of (p,p«) having mean (bt-\, atrit-i), Var[P] = Tf_i, 
Var[n<] =R t = aC t -\ + W tf Corr(p,ji,) « r t = apt- \ Vc M /R<- 



2. Observe the Poisson count y f . 

3. Invoke the Newton procedure of Section 2.1 to obtain estimates of the 
parameters of the approximate posterior distribution of (P,pt) given past 
observations and the new observation 
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To obtain moments for the Poisson mean Xt = fctexp{x/0+pt} note that since 
the distribution of (fit, Pt) is being approximated by a bivariate normal 
distribution, the posterior moments of Xf are approximately 

E[X ( ] * exp xtb t +m t +^[xj x t +C t +2x t p t Jx^Jq\^h t 

Var[X t ] ~ hf exp{2(xtfy + m t ) + [x ( x t +C t + 2x t p t Jxt ]} 

x[exp{^ 2 T ( + C t + 2 x tPt Jx^JC^} * l]. 

The estimate of X<is X t = e btXt+ntt h t . 



2.3 A Simulation Example 

In the example 0 = 0.5, {*/} takes the values {0.25, .5, 1} over and over; ht s i. 

{tod are iid normal with mean 0 and variance 0.25 and a = 0.5. Given 0 and p.*, 

Y t has a Poisson distribution with mean Xt = e Xt ^ + ^ 1 . The simulation starts 

with po drawn from the stationary distribution of {p ( } a normal distribution 

2 1 

with mean 0 and variance W / ( 1-a ) = y Successive p j are computed as 

lit = am-i + ©f; 

the Xt is computed and a Poisson random number with mean Xt is then 
generated. The simulation data are the Poisson counts and {x t }; t = 1, ..., 100. 
The random numbers were generated using LLRANDOM II; cf. Lewis et al. 
(1981). 

At time 0, the Newton procedure is initialized at po = 0, fi = 0, po = h) = 0, 
a = 0.5, W = 0.25, C 0 = 0, r\ = 0.25, To = 1; note the a and W are assumed 
known. The data point y\ is observed and the Newton procedure is used to 
find the posterior moments [m\, b pi, Q, Ti]. 
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The da a point 1/2 is the observed and the Newton procedure is started 
with [m\, b i,r\, C\ + |, Ti ] and used to find [m 2 , £> 2 , pi, C 2 , T2], etc. 

Results of the simulation appear in Figures 1-3. In Figure 1 the count data 
y t appears along with the true X t (dotted line) and the estimated X t (solid line). 
In Figure 2 the true value of }i t (dotted line) and estimated value of fi t (solid 
line) appear. The estimated values of the standard deviation of p t , , also 
appear in Figure 2. The estimated value of /3, P's estimated standard 
deviation yjr t , and the estimated value of p t appear in Figure 3. Figures 1-3 
indicate that the procedure performs satisfactorily. The apparent oscillation 
in some of the figures, (particularly /5 f ), is due to the cyclic nature of {x t }. Not 
surprisingly, the estimated values of X t and p., are less variable than the true 
values. They appear to be practically acceptable. 

3. APPROXIMATE KALMAN PROCEDURES WHICH INCORPORATE 
ESTIMATES OF THE AUTOREGRESSIVE PARAMETERS a AND W: 
NAIVE MOMENT ESTIMATORS 

In this section we will assume that {co t } are independent identically 
distributed normal/Gaussian random variables with mean 0 and variance W 
with (o t independent of {p. ; ; i < f-1), {Y ; ; i< t - 1}, {h t } and p. 

3.1 Estimates of the Autoregressive Parameters a and W 

If [\i t -, t < T } were observable, then one could estimate a and W by 
maximum likelihood; that is, the likelihood function is 




(3.1) 



or 
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(3.2) 



T 

In L(W,a) = l{W,a) = (/i, - a 

2 f =l W 

T T 

— InW — ln27T. 

2 2 



Now 



gives 



and 



from which 



&_ 

da 



T 

-I (ft 



f=l 







= 0 



T T 

a{T) = ]T n t n t .\/ Yj M f 2 -i 



f=l f=l 



a/ 

aw 



1 i 

2 W 2 



£(Mf - «Mf-l ) 2 
f=l 



1^=0 

2 W 




(3.3) 



(3.4) 



(3.5) 



(3.6) 



Unfortunately ^ is not observable. One possible estimate of is the 
posterior mean m t of subsections (2.1) and (2.2). In this case the corresponding 
estimators of a and W are 



“mi?) = J>f-l 

f=l f=i 



(3.7) 






(3.8) 



f-1 
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